# -*- coding: utf-8 -*-

# %%
import matplotlib
import numpy as np
import xlrd
from scipy import ndimage as nd
import time
import os
from matplotlib import pyplot as plt
import scipy.signal
import scipy.io as sio
from scipy.interpolate import UnivariateSpline
import socket

from scipy.optimize import curve_fit

Fit_Prefix = ['NoFit', '']
HiFq_Threshold = 1.1299
Filter_Size = 2
Max_F_Prefix = ['all_freq', 'fq_%.3fEr' % HiFq_Threshold]
Max_F_Switch = 0
Fit_BEC_Switch = 0
Save_Switch = 1

Data_Cut_Shape = 'RectangleRectangle'
font = {'size': 26, 'family': 'Computer Modern Roman'}
matplotlib.rc('font', **font)
# matplotlib.rc('text', usetex=True)

os.system('chcp 65001')
Host_Name = socket.gethostname()
if Host_Name == 'VictorChengPC':
    Computer_Name = 'X:\\Work\\mail.ustc.edu.cn\\Rb87BEC - Files\\'
    os.system('mkdir ' + 'X:\\work\\mail.ustc.edu.cn\\3DSOC_Quench\\figures' +
              Fit_Prefix[Fit_BEC_Switch] + '\\Osc1d\\')
elif Host_Name == 'VictorChengXPS':
    Computer_Name = 'X:\\Work\\mail.ustc.edu.cn\\Rb87BEC - Files\\'
    os.system('mkdir X:\\work\\mail.ustc.edu.cn\\3DSOC_Quench\\figures' + Fit_Prefix[Fit_BEC_Switch] +
              '\\Osc1d\\')
elif Host_Name == 'calcenter2':
    Computer_Name = '\\\\calcenter2\\mail.ustc.edu.cn\\Rb87BEC - Files\\'
# %%
time_start = time.time()
time_start_struct = time.localtime(time_start)


# %%
Lowf_Decay_Expo = 2
Highf_Decay_Expo = 2


def dual_comp(t, t1, a1, f1, phi1, t2, a2, f2, phi2, s0, t0, a0):
    return np.exp(-t**Lowf_Decay_Expo/t1**Lowf_Decay_Expo) * a1*np.cos(2*np.pi*f1*t+phi1) + \
           np.exp(-t**Highf_Decay_Expo/t2**Highf_Decay_Expo) * a2*np.cos(2*np.pi*f2*t+phi2) + \
           np.exp(-t/t0)*a0 + \
           s0


def highf_comp(t, t1, a1, f1, phi1, s0, t0, a0):
    return np.exp(-t**Highf_Decay_Expo/t1**Highf_Decay_Expo) * a1*np.cos(2*np.pi*f1*t+phi1) + s0 + np.exp(-t/t0)*a0


def lowf_comp(t, t1, a1, f1, phi1, s0, t0, a0):
    return np.exp(-t**Lowf_Decay_Expo/t1**Lowf_Decay_Expo) * a1*np.cos(2*np.pi*f1*t+phi1) + s0 + np.exp(-t/t0)*a0


def const_decay(t, t0, a0, s0):
    return np.exp(-t/t0)*a0 + s0


# %%
Data_Table = xlrd.open_workbook(Computer_Name + r'Rb87files\3DSOC_Quench\Quench.xlsx').sheets()[0]
TimeJoin_Table = xlrd.open_workbook(Computer_Name + r'Rb87files\3DSOC_Quench\Quench.xlsx').sheets()[1]
Meta_Table = xlrd.open_workbook(Computer_Name + r'Rb87files\3DSOC_Quench\Quench.xlsx').sheets()[2]
os.system('mkdir D:\\Rb87figures\\3DSOC_Quench\\figuresNoFit\\Osc1d\\')
# MetaRow_List = [19]  # 22860
MetaRow_List = [27]  # 22834

# %%
for MetaRow in MetaRow_List:
    # %% Get data and cal fft
    Cam = Meta_Table.row_values(MetaRow - 1)[5]
    Mat_Data_command = Meta_Table.row_values(MetaRow - 1)[3]
    TimeJoinRow_List = []
    exec(Mat_Data_command)
    Detuning = Meta_Table.row_values(MetaRow - 1)[2]
    if Cam == '22860':
        ky_grid = np.linspace(-1, 1, 66)
    elif Cam == '22834':
        ky_grid = np.linspace(-1, 1, 104)

    Average_SaveFolder = Computer_Name + 'Rb87Data\\3DSOC_Quench\\NoFit\\Selected\\'

    while 1:
        try:
            Data = np.load(
                Average_SaveFolder + 'Cam_%s_Meta_Row=%d,Filter_Size=%s_1d.npz' % (Cam, MetaRow, str(Filter_Size)),
                allow_pickle=True)
            break
        except FileNotFoundError:
            time.sleep(10)

    SP_AllTime_Ave_1d = Data['SP_AllTime_Ave_1d'].item()
    SP_AllTime_Std_1d = Data['SP_AllTime_Std_1d'].item()
    SP_rened_AllTime_Ave_1d = Data['SP_rened_AllTime_Ave_1d'].item()
    SP_rened_AllTime_Std_1d = Data['SP_rened_AllTime_Std_1d'].item()
    ts_TimeJoin = Data['ts_TimeJoin']
    ts = list(SP_AllTime_Ave_1d.keys())
    SP_AllTime_Ave_Array_1d = np.array(list(SP_AllTime_Ave_1d.values()))
    SP_AllTime_Std_Array_1d = np.array(list(SP_AllTime_Std_1d.values()))
    SP_rened_AllTime_Ave_Array_1d = np.array(list(SP_rened_AllTime_Ave_1d.values()))
    SP_rened_AllTime_Std_Array_1d = np.array(list(SP_rened_AllTime_Std_1d.values()))

    SP_FFT_1d = np.fft.fft(SP_rened_AllTime_Ave_Array_1d, axis=0)
    SP_absFFT_1d = np.abs(SP_FFT_1d)
    freqs = np.arange(len(ts)) / (ts[-1] + ts[1]) * 1000
    freqs_Er = freqs / 3.7

    # %% Max peak
    '''Quench_Freq_1d = np.zeros_like(SP_rened_AllTime_Ave_Array_1d[0, ...])
    Quench_Amp_1d = np.zeros_like(SP_rened_AllTime_Ave_Array_1d[0, ...])
    Quench_Freq_HiFq_1d = np.zeros_like(SP_rened_AllTime_Ave_Array_1d[0, ...])
    Quench_Amp_HiFq_1d = np.zeros_like(SP_rened_AllTime_Ave_Array_1d[0, ...])
    y_labels = []
    y_rened = []
    for ky_rened_Index in range(np.shape(SP_rened_AllTime_Ave_Array_1d)[1]):
        ky_Index = ky_rened_Index * Filter_Size
        ky = ky_grid[ky_Index]
        y_rened.append(ky)
        y_labels.append('%.2f' % ky)

        if 0 and 0.3 > ky > -0.3:
            plt.figure(figsize=(16, 8))
            plt.errorbar(ts, SP_rened_AllTime_Ave_Array_1d[:, ky_rened_Index],
                         SP_rened_AllTime_Std_Array_1d[:, ky_rened_Index], fmt='k.-', ecolor='r', ms=15, mfc='c')
            plt.xlim([ts[0], ts[-1]])
            plt.show()

        if Max_F_Switch == 1:
            fq_seled = np.where(freqs_Er < HiFq_Threshold)[0][1:]
        else:
            fq_seled = np.arange(1, len(freqs_Er) // 2 + 1)
        f_ana = freqs_Er[fq_seled]
        SP_absFFT_ana_1d = SP_absFFT_1d[fq_seled, ky_rened_Index]

        Peaks_1d, properties = scipy.signal.find_peaks(SP_absFFT_ana_1d)

        if 0 and 0.3 > ky > -0.3:
            plt.figure(figsize=(16, 8))
            plt.plot(f_ana, SP_absFFT_ana_1d, '-b*')
            plt.plot(f_ana[Peaks_1d], SP_absFFT_ana_1d[Peaks_1d], 'ko')
            plt.plot(f_ana[Peaks_1d[np.argmax(SP_absFFT_ana_1d[Peaks_1d])]], np.max(SP_absFFT_ana_1d[Peaks_1d]), 'ko',
                     ms=20)
            plt.xlim([f_ana[0], f_ana[-1]])
            plt.ylim([0, np.max(SP_absFFT_ana_1d)])
            plt.title('ky=%.3f, Low Freq FFT' % ky)
            plt.show()

        if len(Peaks_1d) != 0:
            Quench_Freq_1d[ky_rened_Index] = f_ana[Peaks_1d[np.argmax(SP_absFFT_ana_1d[Peaks_1d])]]
            Quench_Amp_1d[ky_rened_Index] = np.max(SP_absFFT_ana_1d[Peaks_1d])
        else:
            Quench_Freq_1d[ky_rened_Index] = 0
            Quench_Amp_1d[ky_rened_Index] = 0

        fq_seled_HiFq = np.where(np.logical_and(freqs_Er > HiFq_Threshold, freqs_Er < freqs_Er[-1] / 2))[0]
        f_ana_HiFq = freqs_Er[fq_seled_HiFq]
        SP_absFFT_ana_HiFq_1d = SP_absFFT_1d[fq_seled_HiFq, ky_rened_Index]

        Peaks_HiFq_1d, properties = scipy.signal.find_peaks(SP_absFFT_ana_HiFq_1d)
        if 0 and 0.3 > ky > -0.3:
            plt.figure(figsize=(16, 8))
            plt.plot(f_ana_HiFq, SP_absFFT_ana_HiFq_1d, '-o')
            plt.plot(f_ana_HiFq[Peaks_HiFq_1d], SP_absFFT_ana_HiFq_1d[Peaks_HiFq_1d], 'ko')
            plt.plot(f_ana_HiFq[Peaks_HiFq_1d[np.argmax(SP_absFFT_ana_HiFq_1d[Peaks_HiFq_1d])]],
                     np.max(SP_absFFT_ana_HiFq_1d[Peaks_HiFq_1d]), 'ko', ms=20)
            plt.xlim([f_ana_HiFq[0], f_ana_HiFq[-1]])
            plt.ylim([0, np.max(SP_absFFT_ana_1d)])
            plt.title('ky=%.3f, High Freq FFT' % ky)
            plt.show()

            if len(Peaks_HiFq_1d) != 0:
                Quench_Freq_HiFq_1d[ky_rened_Index] = f_ana_HiFq[
                    Peaks_HiFq_1d[np.argmax(SP_absFFT_ana_HiFq_1d[Peaks_HiFq_1d])]]
                Quench_Amp_HiFq_1d[ky_rened_Index] = np.max(SP_absFFT_ana_HiFq_1d[Peaks_HiFq_1d])
            else:
                Quench_Freq_HiFq_1d[ky_rened_Index] = 0
                Quench_Amp_HiFq_1d[ky_rened_Index] = 0

    fig = plt.figure(1, figsize=(20, 20))
    plt.plot(y_rened, Quench_Freq_1d)
    plt.title('freq')
    plt.show()
    plt.close()

    fig = plt.figure(2, figsize=(20, 20))
    plt.plot(y_rened, Quench_Amp_1d)
    plt.title('Max Peak, Amp')
    plt.show()'''
    # Amp_Ratio_1d = Quench_Amp_1d / Quench_Amp_HiFq_1d

    # %% 2 component fit
    # dual_comp(t, t1, t2, a1, f1, phi1, a2, f2, phi2, s0)
    # sing_comp(t, t1, a1, f1, phi1, s0)
    Freqs_fit = {}
    FreqsLow_fit = []
    FreqsHigh_fit = []
    aLowf_fit = []
    aHighf_fit = []
    FreqsLow_fitStd = []
    FreqsHigh_fitStd = []
    aLowf_fitStd = []
    aHighf_fitStd = []
    y_labels = []
    y_rened = []
    R2_list = []
    Par_list = []
    ParStd_list = []
    ParErr_list = []
    checks_list = []
    for ky_rened_Index in range(26, np.shape(SP_rened_AllTime_Ave_Array_1d)[1]):
        ky_Index = ky_rened_Index * Filter_Size
        ky = ky_grid[ky_Index]
        y_rened.append(ky)
        y_labels.append('%.2f' % ky)
        SP_rened_Ave_Array_1d = SP_rened_AllTime_Ave_Array_1d[:, ky_rened_Index]
        SP_rened_Std_Array_1d = SP_rened_AllTime_Std_Array_1d[:, ky_rened_Index]
        # @@@@@@@@@@@@@@@@@@@@@@@@ Parameters @@@@@@@@@@@@@@@@@@@@@@@@
        # %%
        # @@@@@@@@@@ Single Gauss Parameters @@@@@@@@@@
        if Lowf_Decay_Expo == 2 and Highf_Decay_Expo == 1 and MetaRow == 14:
            Fit_Range = [0, 200]
            p2b2 = np.zeros((11, 3))
            p2b2[0] = [0, 800, np.inf]  # 0, t1
            p2b2[1] = [0, 0.5, 1]  # 1, a1
            p2b2[2] = [0, 1 / 800, 1 / 200]  # 2, f1
            # p2b2[3] = [-np.pi/1, 0, np.pi/1]  # 3, phi1
            p2b2[3] = [-np.inf, 0, np.inf]  # 3, phi1
            p2b2[4] = [0, 500, np.inf]  # 4, t2
            p2b2[5] = [0, 0.5, 1]  # 5, a2
            p2b2[6] = [0, 1 / 100, 1 / 20]  # 6, f2
            # p2b2[7] = [-np.pi, 0, np.pi]  # 7, phi2
            p2b2[7] = [-np.inf, 0, np.inf]  # 7, phi2
            p2b2[8] = [-1, 0, 1]  # 8, s0
            p2b2[9] = [0, 1000, np.inf]  # 9, t0
            p2b2[10] = p2b2[8]  # 10, a0
            if 10 <= ky_Index < 20 or 80 < ky_Index <= 90:
                p2b2[2] = [0, 1 / 800, 1 / 200]  # 2, f1
            if 20 <= ky_Index < 36 or 64 < ky_Index <= 80:
                p2b2[2] = [0, 1 / 1200, 1 / 600]  # 2, f1
            if 36 <= ky_Index < 40 or 60 < ky_Index <= 64:
                p2b2[1] = [0, 0.05, 0.1]  # 1, a1
                p2b2[2] = [0, 1 / 1500, 1 / 1000]  # 2, f1
                p2b2[5] = [0, 0.1, 0.2]  # 5, a2
            if 40 <= ky_Index <= 60:
                p2b2[0] = [0, 2000, 4000]  # 0, t1
                p2b2[1] = [0, 0.03, 0.05]  # 1, a1
                p2b2[2] = [0, 1 / 2000, 1 / 1000]  # 2, f1
                p2b2[5] = [0, 0.1, 0.2]  # 5, a2
                p2b2[6] = [0, 1 / 40, 1 / 20]  # 6, f2
        # @@@@@@@@@@ Double Gauss Parameters @@@@@@@@@@
        if Lowf_Decay_Expo == 2 and Highf_Decay_Expo == 2 and MetaRow == 14:
            Fit_Range = [0, 200]
            p2b2 = np.zeros((11, 3))
            p2b2[0] = [400, 800, 1600]  # 0, t1
            p2b2[1] = [0, 0.2, 0.4]  # 1, a1
            p2b2[2] = [0, 1 / 800, 1 / 200]  # 2, f1
            p2b2[3] = [-np.pi, 0, np.pi]  # 3, phi1
            # p2b2[3] = [-np.inf, 0, np.inf]  # 3, phi1
            p2b2[4] = [50, 200, 800]  # 4, t2
            p2b2[5] = [0, 0.2, 0.4]  # 5, a2
            p2b2[6] = [0, 1 / 100, 1 / 20]  # 6, f2
            p2b2[7] = [-np.pi, 0, np.pi]  # 7, phi2
            # p2b2[7] = [-np.inf, 0, np.inf]  # 7, phi2
            p2b2[8] = [0, 0.5, 1]  # 8, s0
            p2b2[9] = [0, 1000, np.inf]  # 9, t0
            p2b2[10] = p2b2[8]  # 10, a0
            if 10 <= ky_Index < 20 or 80 < ky_Index <= 90:
                p2b2[2] = [0, 1 / 800, 1 / 200]  # 2, f1
            if 20 <= ky_Index < 34 or 66 < ky_Index <= 80:
                p2b2[2] = [0, 1 / 1200, 1 / 600]  # 2, f1
            if 34 <= ky_Index < 38 or 62 < ky_Index <= 66:
                p2b2[0] = [800, 1600, np.inf]  # 0, t1
                p2b2[1] = [0, 0.05, 0.1]  # 1, a1
                p2b2[2] = [0, 1 / 1500, 1 / 1000]  # 2, f1
                p2b2[5] = [0, 0.1, 0.2]  # 5, a2
            if 38 <= ky_Index <= 62:
                p2b2[0] = [800, 1600, np.inf]  # 0, t1
                p2b2[1] = [0, 0.03, 0.05]  # 1, a1
                p2b2[2] = [0, 1 / 2000, 1 / 1000]  # 2, f1
                p2b2[5] = [0, 0.1, 0.2]  # 5, a2
                p2b2[6] = [0, 1 / 40, 1 / 20]  # 6, f2
        # @@@@@@@@@@ Double Gauss Parameters, -0.5Er @@@@@@@@@@
        if Lowf_Decay_Expo == 2 and Highf_Decay_Expo == 2 and MetaRow == 24:
            Fit_Range = [0, 220]
            p2b2 = np.zeros((11, 3))
            p2b2[0] = [400, 800, 1600]  # 0, t1
            p2b2[1] = [0, 0.1, 0.2]  # 1, a1
            p2b2[2] = [0, 1 / 800, 1 / 200]  # 2, f1
            p2b2[3] = [-np.pi, 0, np.pi]  # 3, phi1
            # p2b2[3] = [-np.inf, 0, np.inf]  # 3, phi1
            p2b2[4] = [50, 200, 800]  # 4, t2
            p2b2[5] = [0, 0.2, 0.4]  # 5, a2
            p2b2[6] = [0, 1 / 100, 1 / 20]  # 6, f2
            p2b2[7] = [-np.pi, 0, np.pi]  # 7, phi2
            # p2b2[7] = [-np.inf, 0, np.inf]  # 7, phi2
            p2b2[8] = [0, 0.5, 1]  # 8, s0
            p2b2[9] = [0, 1000, np.inf]  # 9, t0
            p2b2[10] = p2b2[8]  # 10, a0
            if 10 <= ky_Index < 20 or 80 < ky_Index <= 90:
                p2b2[2] = [0, 1 / 800, 1 / 200]  # 2, f1
            if 20 <= ky_Index < 34 or 66 < ky_Index <= 80:
                p2b2[2] = [0, 1 / 1200, 1 / 600]  # 2, f1
            if 34 <= ky_Index < 38 or 62 < ky_Index <= 66:
                p2b2[0] = [800, 1600, np.inf]  # 0, t1
                p2b2[1] = [0, 0.05, 0.1]  # 1, a1
                p2b2[2] = [0, 1 / 1500, 1 / 1000]  # 2, f1
                p2b2[5] = [0, 0.1, 0.2]  # 5, a2
            if 38 <= ky_Index <= 62:
                p2b2[0] = [800, 1600, np.inf]  # 0, t1
                p2b2[1] = [0, 0.03, 0.05]  # 1, a1
                p2b2[2] = [0, 1 / 2000, 1 / 1000]  # 2, f1
                p2b2[5] = [0, 0.1, 0.2]  # 5, a2
                p2b2[6] = [0, 1 / 40, 1 / 20]  # 6, f2
        if Lowf_Decay_Expo == 2 and Highf_Decay_Expo == 2 and MetaRow == 27:
            Fit_Range = [0, 220]
            p2b2 = np.zeros((11, 3))
            p2b2[0] = [400, 800, 1600]  # 0, t1
            p2b2[1] = [0, 0.1, 0.2]  # 1, a1
            p2b2[2] = [0, 1 / 800, 1 / 200]  # 2, f1
            p2b2[3] = [-np.pi, 0, np.pi]  # 3, phi1
            # p2b2[3] = [-np.inf, 0, np.inf]  # 3, phi1
            p2b2[4] = [0, 800, 1600]  # 4, t2
            p2b2[5] = [0, 0.2, 0.4]  # 5, a2
            p2b2[6] = [0, 1 / 100, 1 / 10]  # 6, f2
            p2b2[7] = [-np.pi, 0, np.pi]  # 7, phi2
            # p2b2[7] = [-np.inf, 0, np.inf]  # 7, phi2
            p2b2[8] = [0, 0.5, 1]  # 8, s0
            p2b2[9] = [0, 1000, np.inf]  # 9, t0
            p2b2[10] = p2b2[8]  # 10, a0
            if 10 <= ky_Index < 20 or 80 < ky_Index <= 90:
                p2b2[2] = [0, 1 / 800, 1 / 200]  # 2, f1
            # if 20 <= ky_Index < 34 or 66 < ky_Index <= 80:
            #     p2b2[2] = [0, 1 / 1200, 1 / 600]  # 2, f1
            if 20 <= ky_Index < 32 or 68 < ky_Index <= 80:
                p2b2[0] = [800, 1600, np.inf]  # 0, t1
                p2b2[1] = [0, 0.05, 0.1]  # 1, a1
                p2b2[2] = [0, 1 / 1500, 1 / 1000]  # 2, f1
                p2b2[5] = [0, 0.1, 0.2]  # 5, a2
            if 32 <= ky_Index < 36 or 64 < ky_Index <= 68:
                p2b2[0] = [800, 1600, np.inf]  # 0, t1
                p2b2[1] = [0, 0.03, 0.05]  # 1, a1
                p2b2[2] = [0, 1 / 2000, 1 / 1000]  # 2, f1
                p2b2[5] = [0, 0.1, 0.2]  # 5, a2
                p2b2[6] = [0, 1 / 80, 1 / 30]  # 6, f2
            if 36 <= ky_Index <= 64:
                p2b2[0] = [800, 1600, np.inf]  # 0, t1
                p2b2[1] = [0, 0.08, 0.16]  # 1, a1
                p2b2[2] = [0, 1 / 1500, 1 / 1000]  # 2, f1
                p2b2[5] = [0, 0.1, 0.2]  # 5, a2
        # @@@@@@@@@@@@@@@@@@@@@@@@@@@ Fit @@@@@@@@@@@@@@@@@@@@@@@@@@@
        par, pcov = curve_fit(dual_comp,
                              ts_TimeJoin[Fit_Range[0]:Fit_Range[1]],
                              SP_rened_Ave_Array_1d[Fit_Range[0]:Fit_Range[1]],
                              # sigma=SP_rened_Std_Array_1d[Fit_Range[0]:Fit_Range[1]],
                              p0=p2b2[:, 1], bounds=(p2b2[:, 0], p2b2[:, 2]), maxfev=40000)
        Par_list.append(par)
        x_forfit = ts_TimeJoin[Fit_Range[0]:Fit_Range[1]]
        y_forfit = SP_rened_Ave_Array_1d[Fit_Range[0]:Fit_Range[1]]
        f = dual_comp(x_forfit, *par)
        SS_tot = np.sum((y_forfit - np.mean(y_forfit)) ** 2)
        SS_res = np.sum((y_forfit - f) ** 2)
        R2 = 1 - SS_res / SS_tot
        R2_list.append(R2)
        parstd = np.sqrt(np.diag(pcov))
        ParStd_list.append(parstd)
        ParErr = parstd/np.abs(par)
        if ParErr[2] > 0.2:
            # import pandas as pd
            # df1 = pd.DataFrame(SP_rened_Ave_Array_1d)
            # df1.to_csv("C:\\Users\\cheng\\Desktop\\SP_rened_AllTime_Ave_Array_1d.csv", index=False)
            # df2 = pd.DataFrame(ts_TimeJoin)
            # df2.to_csv("C:\\Users\\cheng\\Desktop\\ts_TimeJoin.csv", index=False)
            time.sleep(0.01)
        ParErr_list.append(ParErr)
        checks = np.zeros((len(par), 3))
        checks[:, 0] = p2b2[:, 0]
        checks[:, 1] = par
        checks[:, 2] = p2b2[:, 2]
        checks_list.append(checks)
        # @@@@@@@@@@@@@@@@@@@@@@@@@@@ Plot @@@@@@@@@@@@@@@@@@@@@@@@@@@
        # %%
        fig = plt.figure(figsize=(25, 15))
        plt.errorbar(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]], SP_rened_Ave_Array_1d[Fit_Range[0]:Fit_Range[1]],
                     SP_rened_Std_Array_1d[Fit_Range[0]:Fit_Range[1]] / np.sqrt(len(TimeJoinRow_List)),
                     fmt='*-', markersize=10)
        plt.plot(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]],
                 dual_comp(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]], *par),
                 'r', linewidth=5, label='dual_comp')

        plt.plot(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]],
                 lowf_comp(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]],
                           par[0], par[1], par[2], par[3], par[8], par[9], par[10]),
                 'g--', linewidth=5, label='low_f')
        plt.plot(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]],
                 highf_comp(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]],
                            par[4], par[5], par[6], par[7], par[8], par[9], par[10]),
                 'm--', linewidth=5, label='high_f')
        plt.plot(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]],
                 const_decay(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]],
                             par[9], par[10], par[8]),
                 'y--', linewidth=5, label='const_decay')

        Start_S = np.ceil(SP_rened_Ave_Array_1d[0] * 10) / 10
        plt.xlim([ts_TimeJoin[Fit_Range[0]], ts_TimeJoin[Fit_Range[1]]])
        # plt.ylim([-0.5, 1])
        plt.ylim([Start_S - 0.6, Start_S + 0.1])
        # plt.title('ky=%.3f,ind=%d,f1=%.3fkHz,f2=%.3fkHz,LfE=%d,HfE=%d' %
        #           (ky, ky_Index, par[2] * 1000, par[6] * 1000, Lowf_Decay_Expo, Highf_Decay_Expo))
        # plt.text(150, 0.9, 'R2=%.3f' % R2)
        plt.text(150, Start_S - 0.1, r'$\tau_1$', fontfamily='Times New Roman')
        # plt.text(150, Start_S - 0.1,
        #          # 'T1=%.2fμs, T2=%.2fμs\n'
        #          r'$\mathbf{f_1}=%.2fkHz*(1±%.3e), \mathbf{f_2}=%.2fkHz*(1±%.3e)$''\n'
        #          'τ1=%.2fμs*(1±%.3e), τ2=%.2fμs*(1±%.3e)\n'
        #          'a1=%.2f*(1±%.3e), a2=%.2f*(1±%.3e)\n'
        #          's0=%.2f*(1±%.3e), t0=%.2f*(1±%.3e), a0=%.2f*(1±%.3e)\n'
        #          'φ1=%.2fπ*(1±%.3e), φ2=%.2fπ*(1±%.3e)\n'
        #          '' %
        #          (
        #           # 1/par[2], 1/par[6],
        #           par[2] * 1000, ParErr[2], par[6] * 1000, ParErr[6],
        #           par[0], ParErr[0], par[4], ParErr[4],
        #           par[1], ParErr[1], par[5], ParErr[5],
        #           par[8], ParErr[8], par[9], ParErr[9], par[10], ParErr[10],
        #           par[3]/np.pi, ParErr[3], par[7]/np.pi, ParErr[7])
        #          )
        plt.text(150, Start_S - 0.55,
                 r'S=$a_1\mathrm{e}^{-t^2/\tau_1^2}\cos(2\pi f_1t+\varphi_1)+'
                 r'a_2\mathrm{e}^{-t^2/\tau_2^2}\cos(2\pi f_2t+\varphi_2)+a_0\mathrm{e}^{-t/\tau_0}+s_0$')
#        plt.legend(loc='upper right')
        plt.xlabel('t/μs')
        plt.ylabel('S')
        plt.tick_params(labelsize=30)
        plt.subplots_adjust(top=0.941, bottom=0.095, left=0.077, right=0.963, hspace=0.2, wspace=0.2)
        plt.show()
        # %%
        if Host_Name == 'calcenter2' and Save_Switch == 1:
            fig.savefig('D:\\Rb87figures\\3DSOC_Quench\\figures' + Fit_Prefix[Fit_BEC_Switch] +
                        '\\Osc1d\\' + 'MetaRow%d,ky=%.3f,ind=%d,f1=%.3fkHz,f2=%.3fkHz,LfE=%d,HfE=%d.pdf' %
                        (MetaRow, ky, ky_Index, par[2] * 1000, par[6] * 1000, Lowf_Decay_Expo, Highf_Decay_Expo))
        elif Host_Name == 'VictorChengPC' and Save_Switch == 1:
            fig.savefig('X:\\work\\mail.ustc.edu.cn\\3DSOC_Quench\\figures' + Fit_Prefix[Fit_BEC_Switch] +
                        '\\Osc1d\\' + 'MetaRow%d,ky=%.3f,ind=%d,f1=%.3fkHz,f2=%.3fkHz,LfE=%d,HfE=%d.pdf' %
                        (MetaRow, ky, ky_Index, par[2] * 1000, par[6] * 1000, Lowf_Decay_Expo, Highf_Decay_Expo))
        elif Host_Name == 'VictorChengXPS' and Save_Switch == 1:
            fig.savefig('X:\\work\\mail.ustc.edu.cn\\3DSOC_Quench\\figures' + Fit_Prefix[Fit_BEC_Switch] +
                        '\\Osc1d\\' + 'MetaRow%d,ky=%.3f,ind=%d,f1=%.3fkHz,f2=%.3fkHz,LfE=%d,HfE=%d.pdf' %
                        (MetaRow, ky, ky_Index, par[2] * 1000, par[6] * 1000, Lowf_Decay_Expo, Highf_Decay_Expo))
        # @@@@@@@@@@@@@@@@@@@@@@@@@@@ Finish @@@@@@@@@@@@@@@@@@@@@@@@@@@
        # %%
        Freqs_fit['%.3f' % ky] = [max(par[2], par[6]),
                                  min(par[2], par[6])]
        if par[2] <= par[6]:
            aLowf_fit.append(par[1])
            aHighf_fit.append(par[5])
            FreqsLow_fit.append(par[2])
            FreqsHigh_fit.append(par[6])
            aLowf_fitStd.append(parstd[1])
            aHighf_fitStd.append(parstd[5])
            FreqsLow_fitStd.append(parstd[2])
            FreqsHigh_fitStd.append(parstd[6])
        else:
            aLowf_fit.append(par[5])
            aHighf_fit.append(par[1])
            FreqsLow_fit.append(par[6])
            FreqsHigh_fit.append(par[2])
            aLowf_fitStd.append(parstd[5])
            aHighf_fitStd.append(parstd[1])
            FreqsLow_fitStd.append(parstd[6])
            FreqsHigh_fitStd.append(parstd[2])
    Par_arr = np.array(Par_list)
    ParErr_arr = np.array(ParErr_list)
    ParStd_arr = np.array(ParStd_list)
    # %% Draw sample(?)
    # fig = plt.figure(figsize=(50, 50))
    # for ky_rened_Index in range(0, np.shape(SP_rened_AllTime_Ave_Array_1d)[1], 5):
    #     plt.subplot(10, 2, ky_rened_Index//5 + 1)
    #     ky_Index = ky_rened_Index * Filter_Size
    #     ky = ky_grid[ky_Index]
    #     y_rened.append(ky)
    #     y_labels.append('%.2f' % ky)
    #     SP_rened_Ave_Array_1d = SP_rened_AllTime_Ave_Array_1d[:, ky_rened_Index]
    #     SP_rened_Std_Array_1d = SP_rened_AllTime_Std_Array_1d[:, ky_rened_Index]
    #     par = Parlist[ky_rened_Index]
    #     ParErr = ParErr_list[ky_rened_Index]
    #     plt.errorbar(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]], SP_rened_Ave_Array_1d[Fit_Range[0]:Fit_Range[1]],
    #                  SP_rened_Std_Array_1d[Fit_Range[0]:Fit_Range[1]] / 2, fmt='*-', markersize=10)
    #     plt.plot(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]],
    #              dual_comp(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]], *par),
    #              'r', linewidth=5, label='dual_comp')
    #
    #     plt.plot(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]],
    #              lowf_comp(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]],
    #                        par[0], par[1], par[2], par[3], par[8], par[9], par[10]),
    #              'g--', linewidth=5, label='comp_1')
    #     plt.plot(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]],
    #              sing_comp(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]],
    #                        par[4], par[5], par[6], par[7], par[8], par[9], par[10]),
    #              'm--', linewidth=5, label='comp_2')
    #     plt.plot(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]],
    #              const_decay(ts_TimeJoin[Fit_Range[0]:Fit_Range[1]],
    #                          par[9], par[10], par[8]),
    #              'y--', linewidth=5, label='const_decay')
    #
    #     Start_S = np.ceil(SP_rened_Ave_Array_1d[0] * 10) / 10
    #     plt.xlim([ts_TimeJoin[Fit_Range[0]], ts_TimeJoin[Fit_Range[1]]])
    #     plt.ylim([Start_S - 0.6, Start_S + 0.1])
    #     plt.title('ky=%.3f,ind=%d,f1=%.3fkHz,f2=%.3fkHz' %
    #               (ky, ky_Index, par[2] * 1000, par[6] * 1000))
    #     #        plt.text(150, 0.9, 'R2=%.3f' % R2)
    #     plt.text(150, Start_S - 0.1,
    #              'T1=%.2fμs, T2=%.2fμs, f1=%.2fkHz*(1±%.3e), f2=%.2fkHz*(1±%.3e)\n'
    #              'τ1=%.2fμs*(1±%.3e), τ2=%.2fμs*(1±%.3e)\n'
    #              'a1=%.2f*(1±%.3e), a2=%.2f*(1±%.3e)\n'
    #              's0=%.2f*(1±%.3e), t0=%.2f*(1±%.3e), a0=%.2f*(1±%.3e)\n'
    #              'φ1=%.2fπ*(1±%.3e), φ2=%.2fπ*(1±%.3e)\n'
    #              '' %
    #              (1 / par[2], 1 / par[6], par[2] * 1000, ParErr[2], par[6] * 1000, ParErr[6],
    #               par[0], ParErr[0], par[4], ParErr[4],
    #               par[1], ParErr[1], par[5], ParErr[5],
    #               par[8], ParErr[8], par[9], ParErr[9], par[10], ParErr[10],
    #               par[3] / np.pi, ParErr[3], par[7] / np.pi, ParErr[7]))
    #     plt.text(150, Start_S - 0.55,
    #              r'S=$a_1\mathrm{e}^{-t^2/\tau_1^2}\cos(2\pi f_1t+\varphi_1)+'
    #              r'a_2\mathrm{e}^{-t/\tau_2}\cos(2\pi f_2t+\varphi_2)+a_0\mathrm{e}^{-t/\tau_0}+s_0$')
    #     plt.legend(loc='upper right')
    #     plt.xlabel('t/μs')
    #     plt.ylabel('S')
    #     plt.tick_params(labelsize=30)
    #     plt.subplots_adjust(top=0.941, bottom=0.095, left=0.077, right=0.963, hspace=0.2, wspace=0.2)
    #%%
    fig = plt.figure(figsize=(20, 20))
    # plt.errorbar(np.array(y_rened), np.array(FreqsLow_fit) * 1000, np.array(FreqsLow_fitStd) * 1000,
    #              fmt='*-k', markersize=10)

    plt.errorbar(np.array(y_rened), np.array(FreqsLow_fit) * 1000,
                 np.where(np.array(FreqsLow_fitStd) / np.array(FreqsLow_fit) < 0.7,
                          np.array(FreqsLow_fitStd), np.array(FreqsLow_fit) * 0.7) * 1000,
                 fmt='*-k', markersize=10)

    plt.plot([-0.3, 0.3], [np.min(np.array(FreqsLow_fit) * 900)] * 2, 'or', markersize=10)
    plt.title('FreqsLow_fit, Detuning=%.3f,LfE=%d,HfE=%d' % (Detuning, Lowf_Decay_Expo, Highf_Decay_Expo))
    plt.xlim([y_rened[0], y_rened[-1]])
    # plt.ylim([0, max(np.array(FreqsLow_fit) * 1000)])
    # plt.ylim([0.5, 1.3])
    plt.xlabel('qy/pi')
    plt.ylabel('f/kHz')
    plt.show()
    if Host_Name == 'calcenter2' and Save_Switch == 1:
        fig.savefig('D:\\Rb87figures\\3DSOC_Quench\\figures' + Fit_Prefix[Fit_BEC_Switch] +
                    '\\Osc1d\\' + 'MetaRow%d,FreqsLow_fit, Detuning=%.3f,LfE=%d,HfE=%d.png' %
                    (MetaRow, Detuning, Lowf_Decay_Expo, Highf_Decay_Expo))
    elif Host_Name == 'VictorChengPC' and Save_Switch == 1:
        fig.savefig('X:\\work\\mail.ustc.edu.cn\\3DSOC_Quench\\figures' + Fit_Prefix[Fit_BEC_Switch] +
                    '\\Osc1d\\' + 'MetaRow%d,FreqsLow_fit, Detuning=%.3f,LfE=%d,HfE=%d.png' %
                    (MetaRow, Detuning, Lowf_Decay_Expo, Highf_Decay_Expo))

    fig = plt.figure(figsize=(20, 20))
    plt.errorbar(np.array(y_rened), np.array(FreqsHigh_fit) * 1000, np.array(FreqsHigh_fitStd) * 1000,
                 fmt='*-k', markersize=10)
    plt.plot([-0.3, 0.3], [np.min(np.array(FreqsHigh_fit) * 900)] * 2, 'or', markersize=10)
    plt.title('FreqsHigh_fit, Detuning=%.3f,LfE=%d,HfE=%d' % (Detuning, Lowf_Decay_Expo, Highf_Decay_Expo))
    plt.xlim([y_rened[0], y_rened[-1]])
#    plt.ylim([0, 8])
    plt.xlabel('qy/pi')
    plt.ylabel('f/kHz')
    plt.show()
    if Host_Name == 'calcenter2' and Save_Switch == 1:
        fig.savefig('D:\\Rb87figures\\3DSOC_Quench\\figures' + Fit_Prefix[Fit_BEC_Switch] +
                    '\\Osc1d\\' + 'MetaRow%d,FreqsHigh_fit, Detuning=%.3f,LfE=%d,HfE=%d.png' %
                    (MetaRow, Detuning, Lowf_Decay_Expo, Highf_Decay_Expo))
    elif Host_Name == 'VictorChengPC' and Save_Switch == 1:
        fig.savefig('X:\\work\\mail.ustc.edu.cn\\3DSOC_Quench\\figures' + Fit_Prefix[Fit_BEC_Switch] +
                    '\\Osc1d\\' + 'MetaRow%d,FreqsHigh_fit, Detuning=%.3f,LfE=%d,HfE=%d.png' %
                    (MetaRow, Detuning, Lowf_Decay_Expo, Highf_Decay_Expo))

    fig = plt.figure(figsize=(20, 20))
    # plt.errorbar(np.array(y_rened), np.array(aLowf_fit), np.array(aLowf_fitStd),
    #              fmt='*-k', markersize=10)

    plt.errorbar(np.array(y_rened), np.array(aLowf_fit),
                 np.where(np.array(aLowf_fitStd) / np.array(aLowf_fit) < 0.7,
                          np.array(aLowf_fitStd), np.array(aLowf_fit) * 0),
                 fmt='*-k', markersize=10)

    plt.plot([-0.3, 0.3], [np.min(np.array(aLowf_fit) * 0.9)] * 2, 'or', markersize=10)
    plt.title('aLowf_fit, Detuning=%.3f,LfE=%d,HfE=%d' % (Detuning, Lowf_Decay_Expo, Highf_Decay_Expo))
    plt.xlim([y_rened[0], y_rened[-1]])
    # plt.ylim([0, 8])
    plt.xlabel('qy/pi')
    plt.ylabel('a1')
    plt.show()
    if Host_Name == 'calcenter2' and Save_Switch == 1:
        fig.savefig('D:\\Rb87figures\\3DSOC_Quench\\figures' + Fit_Prefix[Fit_BEC_Switch] +
                    '\\Osc1d\\' + 'MetaRow%d,aLowf_fit,Detuning=%.3f,LfE=%d,HfE=%d.png' %
                    (MetaRow, Detuning, Lowf_Decay_Expo, Highf_Decay_Expo))
    elif Host_Name == 'VictorChengPC' and Save_Switch == 1:
        fig.savefig('X:\\work\\mail.ustc.edu.cn\\3DSOC_Quench\\figures' + Fit_Prefix[Fit_BEC_Switch] +
                    '\\Osc1d\\' + 'MetaRow%d,aLowf_fit,Detuning=%.3f,LfE=%d,HfE=%d.png' %
                    (MetaRow, Detuning, Lowf_Decay_Expo, Highf_Decay_Expo))

    fig = plt.figure(figsize=(20, 20))
    plt.errorbar(np.array(y_rened), np.array(aHighf_fit), np.array(aHighf_fitStd),
                 fmt='*-k', markersize=10)
    plt.plot([-0.3, 0.3], [np.min(np.array(aHighf_fit) * 0.9)] * 2, 'or', markersize=10)
    plt.title('aHighf_fit, Detuning=%.3f,LfE=%d,HfE=%d' % (Detuning, Lowf_Decay_Expo, Highf_Decay_Expo))
    plt.xlim([y_rened[0], y_rened[-1]])
    # plt.ylim([0, 8])
    plt.xlabel('qy/pi')
    plt.ylabel('a2')
    plt.show()
    if Host_Name == 'calcenter2' and Save_Switch == 1:
        fig.savefig('D:\\Rb87figures\\3DSOC_Quench\\figures' + Fit_Prefix[Fit_BEC_Switch] +
                    '\\Osc1d\\' + 'MetaRow%d,aHighf_fit,Detuning=%.3f,LfE=%d,HfE=%d.png' %
                    (MetaRow, Detuning, Lowf_Decay_Expo, Highf_Decay_Expo))
    elif Host_Name == 'VictorChengPC' and Save_Switch == 1:
        fig.savefig('X:\\work\\mail.ustc.edu.cn\\3DSOC_Quench\\figures' + Fit_Prefix[Fit_BEC_Switch] +
                    '\\Osc1d\\' + 'MetaRow%d,aHighf_fit,Detuning=%.3f,LfE=%d,HfE=%d.png' %
                    (MetaRow, Detuning, Lowf_Decay_Expo, Highf_Decay_Expo))

    # %% First peak
    '''Quench_FirstFreq_1d = np.zeros_like(SP_rened_AllTime_Ave_Array_1d[0, ...])
    Quench_FirstAmp_1d = np.zeros_like(SP_rened_AllTime_Ave_Array_1d[0, ...])
    Quench_FirstFreq_HiFq_1d = np.zeros_like(SP_rened_AllTime_Ave_Array_1d[0, ...])
    Quench_FirstAmp_HiFq_1d = np.zeros_like(SP_rened_AllTime_Ave_Array_1d[0, ...])
    for ky_rened_Index in range(np.shape(SP_rened_AllTime_Ave_Array_1d)[1]):
        ky_Index = ky_rened_Index * Filter_Size
        ky = ky_grid[ky_Index]

        y_rened.append(ky)
        y_labels.append('%.2f' % ky)

        if 0 and 0.3 > ky > -0.3:
            plt.figure(figsize=(16, 8))
            plt.errorbar(ts, SP_rened_AllTime_Ave_Array_1d[:, ky_rened_Index],
                         SP_rened_AllTime_Std_Array_1d[:, ky_rened_Index], fmt='k.-', ecolor='r', ms=15, mfc='c')
            plt.xlim([ts[0], ts[-1]])
            plt.show()

        if Max_F_Switch == 1:
            fq_seled = np.where(freqs_Er < HiFq_Threshold)[0][1:]
        else:
            fq_seled = np.arange(1, len(freqs_Er) // 2 + 1)
        f_ana = freqs_Er[fq_seled]
        SP_absFFT_ana_1d = SP_absFFT_1d[fq_seled, ky_rened_Index]

        Peaks_1d, properties = scipy.signal.find_peaks(SP_absFFT_ana_1d)

        if 0 and 0.3 > ky > -0.3:
            # print(ky)
            plt.figure(figsize=(16, 8))
            plt.plot(f_ana, SP_absFFT_ana_1d, '-b*')
            plt.plot(f_ana[Peaks_1d], SP_absFFT_ana_1d[Peaks_1d], 'ko')
            plt.plot(f_ana[Peaks_1d[0]], SP_absFFT_ana_1d[Peaks_1d[0]], 'ko',
                     ms=20)
            plt.xlim([f_ana[0], f_ana[-1]])
            plt.ylim([0, np.max(SP_absFFT_ana_1d)])
            plt.title('ky=%.3f, Low Freq FFT' % ky)
            plt.show()

        if len(Peaks_1d) != 0:
            Quench_FirstFreq_1d[ky_rened_Index] = f_ana[Peaks_1d[0]]
            Quench_FirstAmp_1d[ky_rened_Index] = SP_absFFT_ana_1d[Peaks_1d[0]]
        else:
            Quench_FirstFreq_1d[ky_rened_Index] = 0
            Quench_FirstAmp_1d[ky_rened_Index] = 0

        fq_seled_HiFq = np.where(np.logical_and(freqs_Er > HiFq_Threshold, freqs_Er < freqs_Er[-1] / 2))[0]
        f_ana_HiFq = freqs_Er[fq_seled_HiFq]
        SP_absFFT_ana_HiFq_1d = SP_absFFT_1d[fq_seled_HiFq, ky_rened_Index]

        Peaks_HiFq_1d, properties = scipy.signal.find_peaks(SP_absFFT_ana_HiFq_1d)
        if 0 and 0.3 > ky > -0.3:
            plt.figure(figsize=(16, 8))
            plt.plot(f_ana_HiFq, SP_absFFT_ana_HiFq_1d, '-o')
            plt.plot(f_ana_HiFq[Peaks_HiFq_1d], SP_absFFT_ana_HiFq_1d[Peaks_HiFq_1d], 'ko')
            plt.plot(f_ana_HiFq[Peaks_HiFq_1d[0]],
                     SP_absFFT_ana_HiFq_1d[Peaks_HiFq_1d[0]], 'ko', ms=20)
            plt.xlim([f_ana_HiFq[0], f_ana_HiFq[-1]])
            plt.ylim([0, np.max(SP_absFFT_ana_1d)])
            plt.title('ky=%.3f, High Freq FFT' % ky)
            plt.show()

            if len(Peaks_HiFq_1d) != 0:
                Quench_FirstFreq_HiFq_1d[ky_rened_Index] = f_ana_HiFq[Peaks_HiFq_1d[0]]
                Quench_FirstAmp_HiFq_1d[ky_rened_Index] = SP_absFFT_ana_HiFq_1d[Peaks_HiFq_1d[0]]
            else:
                Quench_FirstFreq_HiFq_1d[ky_rened_Index] = 0
                Quench_FirstAmp_HiFq_1d[ky_rened_Index] = 0

    fig = plt.figure(1, figsize=(20, 20))
    plt.plot(y_rened, Quench_FirstFreq_1d)
    plt.title('freq')
    plt.show()

    fig = plt.figure(2, figsize=(20, 20))
    plt.plot(y_rened, Quench_FirstAmp_1d)
    plt.title('amp')
    plt.show()'''
    # Amp_Ratio_1d = Quench_FirstAmp_1d / Quench_FirstAmp_HiFq_1d

# %%
time_end = time.time()
time_end_struct = time.localtime(time_end)
print('Time spent: ', time_end - time_start, 's')
print('Start time:', time_start_struct)
print('End time:', time_end_struct)
# sio.savemat('X:\\work\\mail.ustc.edu.cn\\3DSOC_Quench\\figures' + Fit_Prefix[Fit_BEC_Switch] +
#            '\\Osc1d\\' + 'MetaRow%d,SP_rened_AllTime,Detuning=%.3f,LfE=%d,HfE=%d.mat' %
#            (MetaRow, Detuning, Lowf_Decay_Expo, Highf_Decay_Expo),
#            {'Aves': SP_rened_AllTime_Ave_Array_1d, 'Stds': SP_rened_AllTime_Std_Array_1d,
#             'Pars': Par_arr, 'ParStds': ParStd_arr})

sio.savemat(r"C:\Users\cheng\Desktop\quench.mat", {'ts_TimeJoin': ts_TimeJoin,
                                                   'SP_rened_AllTime_Ave_Array_1d': SP_rened_AllTime_Ave_Array_1d,
                                                   'SP_rened_AllTime_Std_Array_1d': SP_rened_AllTime_Std_Array_1d,
                                                   'Par_arr': Par_arr,
                                                   'ParStd_arr': ParStd_arr,
                                                   'ParErr_arr': ParErr_arr})
